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The problem of coupled roll, sway, and yaw stability analysis of submersible 
vehicles is analyzed, with particular emphasis on nonlinear studies. Previous 
results had indicated that a primary loss of stability is through the devel- 
opment of limit cycles. This loss of stability is due to the coupling of roll 
into sway and yaw and cannot be predicted by considering the uncoupled dy- 
namics. In this study, it is shown that the mechanism of loss of stability is 
through bifurcations to periodic solutions. These are characterized as either 
subcritical or supercritical, depending on the sign of a certain nonlinear coeffi- 
cient. Implications of these results to vehicle performance and operations are 


discussed. 
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I. INTRODUCTION 


A. PROBLEM STATEMENT 


The dynamic response of a submersible vehicle operating at the extremes 
of its operational envelope is becoming increasingly important in order to en- 
hance vehicle operations. Traditionally, dynamic stability of motion is studied 
using eigenvalue analysis where the equations of motion are linearized around 
nominal straight line level flight paths [Arentzen & Mandel (1960)], [Clayton 
& Bishop (1982)], [Feldman (1987)]. Directional stability in the horizontal 
plane is normally studied assuming that coupling between sway/yaw and roll 
does not exist. Relaxing this approximation can lead to an oscillatory loss of 
directional stability [Cunningham (1993)] which cannot be predicted by un- 
coupled sway/yaw motions only. This oscillatory loss of stability can generate 
limit cycles in the system, as was confirmed numerically in previous studies 
[Cunningham (1993)]. In order to gain a better understanding of the mecha- 
nism of this type of loss of stability and the stability properties of the resulting 
limit cycles, it is necessary to perform a systematic nonlinear analysis which 


is precisely the scope of this work. 


B. OBJECTIVES AND OUTLINE 


In this work we examine the problem of stability of motion with controls 
fixed in the horizontal plane, with particular emphasis on the mechanism of 
loss of stability of straight line motion. Coupling between sway/yaw and roll 
is taken into consideration. This has its origins in both inertial and hydro- 
dynamic coupling. We concentrate on an oscillatory loss of stability case, 
where one pair of complex conjugate of eigenvalues of the linearized system 
matrix crosses the imaginary axis. This loss of stability occurs in the form 
of generic bifurcations to periodic solutions [Guckenheimer & Holmes (1983)]. 
Taylor expansions and center manifold approximations are employed in or- 
der to isolate the main nonlinear terms that influence system response after 
the initial loss of stability [Hassard & Wan (1978)]. Integral averaging is 
performed in order to combine the nonlinear terms into a design stability co- 
efficient [Chow & Mallet—Paret (1977)]. Special attention is paid to the study 
of the quadratic drag terms as they constitute some of the main nonlinear 
terms of the equations of motion. The difficulty associated with the nons- 
moothness of the absolute value nonlinearities is dealt with by employing the 
concept of generalized gradient [Clarke (1983)], a technique which was utilized 
in [Papadimitriou (1994)]. This has the advantage of keeping the linear terms 
constant, unlike the linear /cubic approximation typically used in ship roll mo- 
tion studies [Dalzell (1978)], where the linear damping coefficient is a function 


of the assumed amplitude of motion. 


Vehicle modeling in this work follows standard notation [Gertler & Hagen 


(1976)], [Smith et al (1978)], and numerical results are presented for a variant 
of the Swimmer Delivery Model used in [Cunningham (1993)] for which a set 
of hydrodynamic coefficients and geometric properties is available. Although 
the main results and conclusions of this work are derived for a submerged 


vehicle, similar techniques can be applied to surface ships as well. 





II. EQUATIONS OF MOTION 


A. COORDINATE SYSTEM 


In our analysis we are going to use a moving coordinate system (z, y, z), 
attached on the vehicle. The origin of this system coincides with the center 
of buoyancy, B. The z-axis is attached to the longitudinal plane of symmetry 
for the vehicle, the y—axis is positive starboard, and the z—axis is positive 
downwards. All main symbols used in the development of the equations of 


motion in this work are summarized in ‘able 1. 


B. GENERAL FORM OF THE EQUATIONS OF MOTION 


The equations of motion for a submersible vehicle in the horizontal plane 


are written as follows: 


Sway equation: 


mo + Ur — wp + z¢(pq + 7) — yo(p* +r) + 2e(ar — )] = 
DGD Glee DCist Sg Gime te Oe yOt 

Y5U75, + Yod + VYpUp + Y,Ur + Yogug + Yupwp + Yurwr + 
(W — B)cos@sing — 


[™ [eo,n(2)o + 27)? + Cp.v(@)(w - 29)" ane) 


tail Uo oe () 


Yaw equation: 


Iz? + (yy — Inz)pq — Iey(p* — 9°) — Iyz(pr + 4) + Inz(qr — p) + 
m|ze(o + Ur — wp) - yg(U —vr+wq)] = 
Nap + Nzr + Nogpg + Ngrar + N.Uv + Nowow + 
N5,U76, + Ng + NpUUp + N,-Ur + Nogvg + Nupwp + Nupwr + 
(cqgW — rpB)cos@sing + (ygW — ypB) sin + U*Nprop — 

V+ ie 
are 


i [Cp,h(x)(v + ar)? + Cp,b(x)(w — x4)"| ( 


tail Daal ie) oa (2) 


Roll equation: 


InaB + (Ize — Iyy)ar + Iny(pr — 4) — Iyz(q? — 7”) — Iez(pq +7) + 
mlyg(w —Uq+ vp) — zg(v + Ur — wp)| = 

Kap + Ker + Kyopg + Korgr + A,Uv + Kyyvw + 

K propU* + Kg + KpUp + K,Ur + Kygvg + Kupwp + Kypwr + 


(ygW — ypB)cosécos¢ — (zgW — zgB)cos@sing . (3) 


The rotational velocity equation around the x-axis is simply, 
P=D, (4) 
while U.¢ denotes the cross-flow velocity, 


Ucs(x) = yi (v + ar)? + (w — 29). 








translational velocities about (z, y,z) axes 





=e es es es ee Gee es) ee ee ee es Gees Gey 





Table 1: Nomenclature 
C. SIMPLIFICATIONS 


Before we proceed with the analysis, we must simplify the above equations 
of motion in order to reflect the fact that we are analyzing motion in the 


horizontal plane only. The simplifications that we employ are the following: 


Velocity, w, and acceleration, w, in the z-direction are zero. 


e Acceleration in the longitudinal direction, wu, is zero. 
e Rotational velocity, g, and acceleration, g, in the y-direction are zero. 


e Pitch angle, @, is zero. 


D. SIMPLIFIED EQUATIONS OF MOTION 


After the application of the above simplifying assumptions, the equations 


of motion become: 
Sway equation: 


(m — Y,)v + (mag — Y;)f — (mzg + Yp)p = 
Y,Uv + (Y¥,U — mU)r + YpUp + yep" + yer* + (W — B)sing + ¥s,U75, — 


? Cp,h(2)(v + aay (vu +czr) 


. Tae (5) 


Yaw equation: 


(mzG a Ng)v ar ls = N;)r = (las a N5)p = 
N Uv + (N,U — mzgU)r + NpUp + Inyp* + Iyzpr + your + 
(rgGW _ rpBB) sin @ + N5,U*6, A a poe _ 


Lnose (v + rr) 
Cp. h(z\@ + cr)? ——__— 
i. ae ( ) ele) 


a9 ge | (6) 
Roll equation: 


(—mzg — Ky)v + Uz — Ke)? + Ur — Ks) p = 
K,Uv + (K,U + mzgU)r + KpUp — Ipypr — Tyr? — myqup + 


(ycW — ypB) cos¢ — (zgW — zgB)sind+ U*Kprop - (7) 


Roll rate: 


CO 


Ill. LINEAR ANALYSIS ~. 


A. LINEARIZATION 


The simplified equations of motion can be written in matrix form as follows: 
Ax =Bx+g(x), (9) 


where the state vector x is defined as, 


OSs se] 


and the state matrices are, 


m—Yz HIME GE VES ANGE — NG 
mice — Nagle, Ne ie, — 1X 
=e ig ies a kee / (ee oe K; 
0 0 0 


= 


— ==] © 


and 

YGo YU — mu Y,U 

N,U -—mzgGU+N,U N,v 

KU Siecle he 
0 0 1 


ie = 


oS oS 


The term g(x) contains all the nonlinear terms, 


g1 = yop? +ycr? +(W — B)sing + ¥s,U6, — 


Cp, / om h(x)(v+ar)|u+c2r|dz , (10) 


tail 


| 


92 Ieyp” eed GUT (rgW a rpB) sin @ + N¢,U76, = 


Cp, [ h(xz)(u+ar)|v + cr|z dz , (11) 


G35 ae —Ipypr = et — WMG Eup 1 USS. a 


z 


(yGW — ypB)cos¢ — (zgW —zpB)sing, - (12) 


CA QO. (13) 


We want to linearize the nonlinear terms about a nominal point 


Xo = [v9, 70, Po, do)’ = 0. 


After applying Taylor series expansion of the non-linear terms about the nom- 
inal point xp and keeping only the linear components, the linearized equations 


of motion are written in matrix form as: 
A'x =Bx, (14) 


where 
Ee 


and 
au Y,U —mU ie W-B 


pee mocl +0V,0 NiU “gew —agB 
Poyometiceu 11,0 8i,U  —726W + Zeb 
0 0 1 0 


Equation (14) is our linearized system. Eigenvalue analysis for this system is 


= 


performed in the next section in order to assess the dynamic stability of the 


vehicle. 


B. LOSS OF STABILITY 


Stability of the linearized system depends on the location of the four eigen- 


values of the system, in other words the roots of the characteristic equation: 
dete —— A 0), ells) 


10 


or, in polynomial form, 


AM+ Bi +CN+4+ DA+E=0. (16) 


The coefficients of the polynomial equation (16) are given by: 


A 


(Inc — Kp)[(m — Ys)(Izz — Nz) — (mze@ — Yz)(mrg — N»)] 

(mzg + K,)|(—I. + N;)(mze + Y,) + (NV; + 1..)(mae — Y;)] 
(A;)[(~mag + Ns)(mzg + Ys) + (Np + Iez)(m — Yo) , (17) 
KU|(—Iez + Ne)(mzg + Ys) + (Np + Inz)(mag — Y3)] 

(mzg + Ky)[(zz — N+)(YpU) — (UN, — Umag)(mzg + Y5) 

(Np + Iez)(UY, —- Um) — (NpU) (mae — Y;)] 

(Ina — Ky)[(m — Ys)(UN, + Umag) + (YU) (Lez — N3) 

(mag — Y:)(N.U) — (mag — Ny)(UY, — Um)] 

KpU[(m — Ys)(Izz — Nr) — (mag — Y;)(mag — Ny) 

(K; — Inz)|(mag — Ng)(YpU) — (NyU)(mzg + Yp) 

(Np + Izz)(YoU) — (NpU)(m — Ys) 

(mzgU + K,U)|(—mag + Ns)(mze + ¥5) 

(Np + Izz)(m — Y5)] (18) 
(Ine — Ky)((YoU)(UN, — mUzg) — (UY, — Um)(N.V)] 

K,U[(m — Y5)(UN, — mUzg) + (YU) (Lez — Ni) 

(mag — ¥:)(NyV) — (mag — Ns)(UY; — Um) 

(mzgU + K,U)|(mag — Nu)(YpU) — (NoU)(mze + Yop) 

(Np + Inz)(YoU) — (NpU)(m — Y5)] 


i 


+ (zgW — 2pB)|(m — Y5)(Izz — Ns) — (mag — ¥s)(mag — Ns) 

— (mzgo+t+ K;)|(ceB — rtgW)(marg —- Y;) - (UN, — Umzg)(YpU ) 

+ (UY, —Um)(N,U)) 

— K,U[Uz2 — N+)(YpU) — (UN, — Umzg)(mzg + Y5) 

+ (Nz + Izz)(UY, — Um) —(NpU)(mze — Y3)] 

+ (K;—Izz)[(teaB —tgW)(m — Y;) 

— (NU)(¥,U) + (NU)(YU)) (19) 
D = (mzgU + K,U)|(r4pB — cgW)(m — Y5) 

— (NyU)(¥pU) + (NpU)(YU)) 

~ (K;—Is:)(zaB - zgW)(YU) 

+ (mzg+Ky)(zgB —xgW)(UY, — Um) 

— K,U|(cgB - cgW)(mzg — Y;) 

— (UN, —mUzg)(YpV) + (UY, — Um)(N,U)] 

~ KpU[(YU)(UN, — mU ag) - (UY; — Um)(NwU)] 

— (zgW - zpB)|(m — ¥;)(UN, — mUzg) + (YoU)(Iee - Ne) 

— (mzg — ¥;)(NvU) — (mag — Ns)(UY, - Um)), (20) 
E = (zqW —- zgB)[(YU)(UN, — mUzc) — (UY, — Um)(N,U)] 

+ (K,U)\(2pB — zgW)(UY, — Um) 


— (mzgU + K,U)(x#pB — rtgW)(Y,U) (is 


We can examine the stability of the system by utilizing Routh’s criterion. 


Application of this criterion to the characteristic equation (16), reveals that 


TZ 


the following two conditions must be satisfied in order to ensure that all roots 


of (16) have negative real parts: 


BOD Ay ane? =). (22) 


E>OQ. OB) 


If E is less than zero, one real root of (16) becomes positive and the sytem will 
become unstable in a divergent manner (Guckenheimer and Holmes, 1983). 
This is the case of a directionally unstable ship which is well known in the 
literature (Clayton and Bishop, 1982). If, however, condition (22) is violated, 
the system will exhibit an oscillatory motion due to the presence of complex 
conjugate roots with positive real parts. This form of instability is caused by 


the coupling of roll into sway and yaw and is further analyzed in this work. 


In order to compute the limiting case of loss of stability, we consider equa- 
tion (24), 
BCD — AD* —- EB*=0. (24) 


The result of this equation will produce the limiting value of zg as a function 


of xg for loss of stability. A curve of this functional form, 


zg = f(z@) , 


will be our locus of loss of stability. After some algebra, we can express the 


coefficients of equation (24) in the following form: 


A= Aize + Aozgot+ Az, (25) 


13 


where 


ro 


Ag = 


A3 = 


where 


oi 


—m?*(Izz — N;) 

—mY,(I22 — Nz) — mMK5(Izz — Nz) + m(Nz + Inz)(maeg — Y;) 
m(K; — Izz)(mzg — Nz) 

(Ine — K5)(m — Ys)(zz — Ne) — za — Kp)(mag — Y;)(mze — Nz) 
Ig all IN) EG GEE Cipere: — 0) 


ne = Pe eves = N3) or CG ic Tx2)(No5 a5 Iz2)(m _ Ys) 


B= Bizz + Boze + Bs, (26) 


m*(UN, — Umzg) + m*U(mzg — Ny») 

—m( Ky U)(Iz2 — Nz) — mM(Uz2 — Nz)(YpU) 

mY,;(UN, — Umzg) + mK,(UN, — Umzg) 

m(Np + Izz)(UY; - Um) + m(N,U)(mzg — Y;3) 

m(K; — Izz)(N,U) + mUY;(mag — Ny) 

mU (N53 + Inz)(m — Ys) + mUK,(mzg — Ny) 

~Y;(KyU)(Iez — Ns) + (KU)(Np + Inz)(mag — Ys) 

Ky(Izz — Nz) (YpU) + KsY,(UN, — Umze) 

Ky(Ng + Inz)(UY; - Um) + Ki(N,U)(mze — Y;) 

(Ure — Kp)(m — Yy)(UN, — Umag) — Ure — Kp)(YoU) (Tez — Ne) 


(Ine — Ky)(mag — Y;)(NyU) + Use — K5)(mag — Ny)(UY, — Um) 


14 


where 


C7 


Co 


C3 


(K,U)(m — Yi)(Is2 — Nz) + (KpU)(mae — ¥s)(mag — Ne) 
(K; — Ine)(mag — No)(¥U) ~ ¥(Ks — Inz)(NoU) 

(Np + Iaz)(YsU)(Ky ~ Ine) — (Ki — lez)(NpU)(m = Ys) 
UK,Y;(mze@ — No) 


C= Gis + Cozet C3, (27) 


—m°U(N,U) 

mU (mzg — Ny)(YpU) — mUY;(N,U) — mUK,(N,U) 

mU (Np + Inz)(YyU) — mU(NpU)(m — Y;) 

W (m — Ys)(Izz — Nz) — W(mag — Yz)(mag — N3) 

m(XpB —atgW)(mzg — Y;) + m(UN, — Umzg)(¥YpU) 
m(UY, — Um)(N,U) + m(K,U)(UN, — Umzg) 

(Ing — K3)(YoU (UN, — Umag) — (Ine — Kp)(UY; — Um)(N,U) 
(K,U)(m — Y¥s)(UN, — Umzg) + (KU) (YU) (Lez — Nz) 
(KpU)(mae@ — Yz)(NU) — (KpU)(mag — Ny)(UY, — Um) 
UK,(mrg — Ns)(YpU) — UK;Y5(NwU) 

U K,(Ng + Inz)(YoU) — UK,(NpU)(m — Y5) 

Ki(XpB —xzgW)(mzg — ¥;) + Kx(UN, — Umag)(YpU) 
K(UY, — Um)(NpU) — (Ky) (Izz — Ne)(YpU) 
Y;(KyU)(UN, — Umag) — (KyU) (Np + Inz)(UY, — Um) 
(KyU)(NpU )(mag — Y;) + (Ki — Inz)(XBB -— teW )(m — Yo) 


1s) 


Se — a) (NU YG ) Ke —Isz)( Np) 
De= ize Ds. (28) 
where 


Di, = mU(xgB —xgW)(m — Y;) — mU(N,U)(YpU) 
+ mU(N,U)(Y,U) + m(xpB — zgW)(UY, — Um) 
~ W(m-—Y;)(UN, — Umzg) — W (YU) (Ize — Ne) 
+ W(mzg — Y;)(NyU) + W(mzg — Ns)(UY; — Um) 

Dy = UK,(zpB- rgW)(m — ¥;) —UK,(N,U)(YpU) 
+ UK,(NpU)(Y.U) — (Kp — Iez)(agB — tgW)(YU) 
+ Ky(cgB —xgW)(UY, — Um) — (K,U)(zgB — tgW)(mzg — Y;) 
+ (K,U)(UN, — Umag)(¥pU) — (K,U)(UY, — Um)(NpU) 
— (K,U)(Y,U)(UN, — Umzg) + (KpU)(UK, — Um)(N,U) 

and 
E = E,zg+ Eo, (29) 


where 


FE, W(Y,U)(UN, —Umzcg) — W(UY, — Um)(N,U) 
— mU(rpB-xrgW)(Y,VU) 


Ey = (K,U)(xpB —zrgW)(UY, —Um) —UK,(rzgB —xgW)(Y,VU ) 


If we apply the stability criterion (24) utilizing expressions (25) through 


(29), we get a fifth order polynomial equation in the metacentric height z¢ of 


16 
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Figure 1: Critical value of zg versus xg for U = 5 ft/sec 
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Figure 2: Critical value of zg versus xq for different values of U (ft/sec) 
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the following form, 
Fsza t+ Faza + F3z3 t+ Foz + Fizg + Fo = 0, (30) 
where, 


Fy = B3C3Do — A3D3 — E2B?, 

F, = B3C3D, + (B3C2 + BoC3)D2 — 2E2B2B3 — E\B? - 
2D2D1A3 — D3A2 , 

Fp = —E(B} + 2B,B3) — 2E,B2B3 + (B3C2 + BoC3)D, + 
(B3C, + BoC2 + BC3)D2 — D?.A3 — 2D2D1 Aq — D2 A, , 

F3; = —D?A,—2D2D,A, — 2E2B,B2 — E:(B2 + 2B,B3) + 
D1(B3C; + BoC. + BiC3) + Do(BoC; + BiC2) , 

Fy = D,(BoC, + ByC2) + B,C, D2 — D?A, — E,B? — 2E,B,Bz , 


F; = —E,B?+B,C,D}. 


Using the corresponding Matlab program shown in Appendix A, we can solve 
equation (30) using a typical constant forward velocity U = 5 ft/sec. The value 
of xg is given in the range from 0 to 2 ft. In this way we get five solutions in 
zg for a given value of xg. Investigating these solutions, we can see that only 
one satisfies the second criterion (23). This solution corresponds to the locus 
of loss of stability. For values of zg greater than the critical value, the system 
is stable, and for values less than its critical value, the system is unstable, see 
Figure 1. The calculations are repeated for a range of forward speeds U from 


2 to 8 ft/sec. The results are shown parametrically in Figure 2. We observe 
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an upwards movement of the curves as the speed is increased. This means 
that the system experiences a tendency to become less stable at higher speeds 


which, in turn, calls for higher metacentric heights to ensure stability. 


1 





IV. NONLINEAR ANALYSIS 


A. INTRODUCTION 


From the linearized analysis on the loss of stability that was done in the 
previous chapter, we can see that as a specific parameter of the system, such as 
(xg, zg), is varied, it is possible to pass from a region of stability to a region of 
instability. This case of loss of stability is associated with one pair of complex 
conjugate eigenvalues of the system crossing the imaginary axis. This loss of 
stability is usually accompanied by self sustained oscillations, and it is called 


Hopf Bifurcation. 


There are two cases of Hopf Bifurcations, supercritical and subcritical. In 
the supercritical case, limit cycles are created just after the loss of stability. 
A limit cycle is a constant amplitude oscillatory motion of our system. This 
amplitude is usually larger as we move away from the bifurcation point. When 
the limit cycles have small amplitudes the situation is not very critical for our 
vehicle and is not much different from stable conditions. In the subcritical 
case we may have convergence to limit cycles, even before our system looses 


its stability. Furthermore, the limit cycle amplitudes are considerably higher. 


The analysis that will follow will be performed in order to verify the exis- 
tence of the limit cycles and to find out which case of Hopf Bifurcations we 
have in our model. We will also examine the stability of the resulting limit 


cycles. This analysis is necessary because we can predict the behavior of a 
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vehicle when some of its parameters change and we have operation in the 
proximity of a bifurcation point. The results of the non linear analysis will 
be verified by a numerical simulation of the system’s motion, both above and 


below the bifurcation point. 


B. THIRD ORDER EXPANSIONS 


From the linearization procedure we see that our system is written in the 
form of matrix equation (14), where we have ignored the non linear terms. 


If we take into account the non linear terms up to third order, equation (14) 


becomes: 
A’x =B’x+ g(x) ; (31) 
where, 
ane 
_ | ga(x) 
g(x) 93(x) ’ 
ga(x) 


Keeping terms up to third order the vector of non linear terms can be written: 
g(x) = g(x) + g(x) + (x) , (32) 


where g'?)(x) contains the second order non linear terms, g)(z) contains the 
third order non linear terms, and c(xr) contains the constant terms. The cross 


flow integrals can be written as follows: 


‘ey Co, | ~ h(x)(v + a2r)|v + zr| dz 


tail 


Cp, | _ h({x)(u+ar)|v + or|x dz 


tail 


a 
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(33) 


The second order non linear terms are: 


a2) 
(2) a 9 (8) 
& (x) = 2 ? 
05. iG 
94 (x) 
where 
Q = yep? + yor? — 1) 
on = Inyp* + Iyzpr + your — 1 
(2) iad ay : 2am _ i 2 
93 — xzyPT yzl MYGUr (ygW ypB )o 
2 
9° = (0. 
The third order non linear terms are: 
3 
2. (x) 
4 6 
g(x) = 2 <i 
a (Ge) 
93 (z) 
where, 
1 
Qe 5 iW — B)¢° 
1 
gg = —-18) — a(toW — 2pB)¢° 


1 
gs = 5 (2eW — 2pB)¢° 


3 
9 = (0. 


The constant terms are: 


where, 
co(z) = Ng,U*5, + U*Norop 


SS 
oO 
—-_~ 

8 
ee 

| 


c4(x) = Oe 


In order to get the second and third order non linear terms of the cross flow 


integrals, we must expand in Taylor series a function of the form: 


‘flO =éle 


about a nominal point €o: 


EE] = Eoléo| + 2lé0l(€ — £0) + sign(Eo)(€ — £0)? + fP(E) . 


The sign function in equation (33) can be approximated by: 


sign(o) = lim tanh(22) , 


where the approximation gets better as 7 gets smaller. 


If we choose 9 = 0 as our nominal point, equation (35) becomes: 


l 43 
é\éi = 6° 
In our case € is (v + zr) so we have: 
1 3 
(u+ear)|v+ar| = —(v4+cr) 
Oy 


OT 


it 
(u+ar)|v+azr| = mle + 2°r? + 30% cr + 327r7v) 
a 
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(34) 


(35) 


(37) 


(38) 


(39) 


Using equation (39) the cross flow integrals become: 


where, 


D 





@ 
I, = ¥(Egv® + 3E,v2r + 3Equr? + E3r°) 


Oy 


C 
ie = 5 (E10 =F 3Eov-r Ie 3E3ur" aE E,4r°) 
y 


1p -_ t'*h(x)dz ,i =0,1,2,3,4 


Ltail 


(40) 


(41) 


(42) 


By using this approximation we see that the expansion of the cross flow inte- 


grals give only third order terms, so we have: 


12) — J) = 9 


(43) 


Now we want to write our matrix equation in state space form, so we have to 


find the inverse of the system matrix: 


where, 


Q11 
Q12 
43 
Q21 
22 
23 


Q31 


1 42 213 Q 
a le?) 
O21 222 923 Q 
(Ayt=| 2) 2 P 
@31 a32 33 0 ) 
i YD 
0 Oa 


(Iez — N+)(ea — Kp) — Uz — Ke)(—Inz — Np) 

(Y¥;— mzq)(Uze — Kp) + zz — K+)(—mze — Y5) 
(mag — Yr)(—Izz — Np) — Uz22 — Nr)(—mze — Y5) 
—(mzg — Ns) (Tez — Kp) + (—mze — Ko)(—Iez — Ny) 
(m — Ys) (In — Ky) — (—mzq — Ks)(—mze — Y5) 
—(m — Ys)(—Inz — Np) + (mag — No)(—mze — Ys) 
(mag — Ny) (Izz — K+) — (—mzg — Ka) (Tez — Ne) 
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a3q = —(m—Y%)U2z — K;) + (—mzg — Ky)(mzeg — Y;) 
233 >= (m — AOE = N;) — (mz¢ — Ny)(mze¢ ae a) 
and 
( (m = ayes oy N;) Lee a K%5) 
eee) (=i — IN) 
= (mz¢ ae Ns)(mz¢ — Veer = le) 
+ (mzg — Ny)(Isz — K;)(—mze — Y5) 
== (—mzc¢ bez K,)(mze¢ = Y;)(—Izz a N35) 
— (-—mzg — Ky)(Izz — Nz)(—mzeg — Y5) 
If we multiply equation (31) by (A’)7} from the left side we get: 
(A) A'g = (A’)7* B's + (A')7*9(z) 
or 
1p Wi ie 1G fas) 


where, 


LNG So, ee ae Cl 
Fn Fo Fax Fas 
P=) A Fo fa Fu | > 
joe 0) oe 
Cm 0 4] we 


(44) 


(46) 


(47) 


where, 


az(Y,U ) + ajo(NU) + ai3( KV) 

aii(Y, — m)U + ayo(N, — mag)U + a33(K; + mzg)U 
ayi(YpU) + ai2(NpU) + a13(K pV) 

aui(W — B) +a 2(rxgW — rpB) + 413(—zgW + zBB) 
aa (YyU) + ao2(NyU) + ao3( KU) 

aai(Y, — m)U + age(N; — mzg)U + ao3(K, + mzg)U 
a21(YpU) + a22(NpU) + a23(KpU) 

an (W — B) + a22(zGW — xpBB) + a23(—zGW + zpB) 
a3i(Y,U) + a3o( NU) + a33(K.U ) 

a3i1(Y¥; — m)U + a3o(N, — mzg)U + a33(K, + mzg)U 
a31(YpU) + a3e(NpU ) + a33(KpU) 


a3i(W — B) + a3o(rxgW — xBB) + 433(—zGW + zBB) 


From equations (45),(47) we see that: 


where, 


Gi(z) 
G(x) =(4)o(a) = | B27 | 
Ga(z) 
Giz) = Bole) + S9alz) + F99(2) 
Gz) = Fare) + Foal) + F9s(c) 
Go(z) = Sor(z) + oa(2) + F-9a(2) 
Gi(z) = 0 


Zit 


Equation (45) can also be written in the following way: 
e=Fre+G(r)+G°(c) +k (48) 


where, 


GA (x 


and 


Each element of the above non linear terms vectors can be written as follows: 


2 Qi1 (2 Q12 (2 CTS (2 
GP'(e) = gi(2) + 93 (2) + $95" (2) 
2 Q221 (2 Q22 (2 223 (2 
GP (2) = Raia) + of (c) + Bac) 
2 Q31 (2 Qa32 (2 Q33 (2 
Gs (x) = 5 (zx) aE 92 (x) a ~ 93 (zx) 

G?)(z) = 0 
and 
3 Qi1 (3 Q1i2 (3 Q13 (3 
GPa) = ar (2) + 93) (2) + 9$"(e) 
3 Q21 (3 Q22 (3 Q23 (3 
E2(@) = 9 (x) + 93 (x) + 93 (x) 
3 Q31 (3 A232 (3 a33 (3 
G})(z) = er (x) + 92 (x) + = 93 (x) 
G2(z) = 0 
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Finally, the constant terms are: 


Qj1 Q12 Q13 


os eee GoD 
Q21 a22 223 
| me= exe 
2 D C1 + D C2 ae D C3 
Q31 a32 233 
SO le pea + Dp? + DS 
See oi ae UB 


C. COORDINATE TRANSFORMATIONS 


From equations (45) and (48) it is obvious that the stability of our system 
depends on the eigenvalues of matrix F. Since we want to investigate the 
behavior of our system oround the Hopf Bifurcation point, it is useful to bring 
our system in its normal coordinate form. This can be done by applying a 
transformation of the coordinate system, using as transformation matrix T, 
the modal matrix of eigenvectors of F’, evaluated at a critical point. This 
critical point will be a pair of zg and zg values, that belong to the critical line 


of loss of stability. The applied transformation will be as follows: 
ip agi (49) 


Or, 
z=T 2 (50) 

Then equation (48) can be written as follows: 
Tz = FTz2+G (Tz) + G@(Tz) +k (51) 


DNS, 


Or, 


= Te Pie PAG (Tz) metG (Pz) +7 tk (52) 


In this new coordinate system, the system’s matrix is T~’ FT, and at the Hopf 


Bifurcation point can be written as follows: 


0 =a Oe 
Pe ltie 0 Oe 

ee 0 eae 
0 0 OP 


Where wp is the imaginary part of the critical pair of eigenvalues, and the 
remaining eigenvalues P; , Pp. are negative. For values close to the Hopf Bi- 


furcation line, we can write the system matrix as follows: 


a’e —(wo + w’e) 0 0 
tng _ | (wot+w’e) a’e 0 0 
= 0 0 (Py+Ple) 0 : 
r 0 dD | esse) 


where: « = Difference from the critical point (zg — z¢.). 
a’ = Derivative of the real part of the critical eigenvalue with respect to e. 


w’ = Derivative of the imaginary part of the critical eigenvaluewith respect 


to €. 
P| = Derivative of P; with respect to e. 


P3 = Derivative of Pz with respect to e. 


D. CENTER MANIFOLD EXPANSIONS 


By using the coordinate transformation described in the previous chapter, 
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and from equation (50), we see that we have a new variables vector, which is: 

z] 

zo 

23 

24 
The first two coordinates z, , z2, are the critical coordinates and correspond 
to a pair of complex conjugate eigenvalues. The remaining two coordinates 
z3, 24, are the stable coordinates and they always correspond to eigenvalues 
that are negative. The center manifold theory predicts that the relationship 
between the critical coordinates z), z2, and the stable coordinates z3, z4, is at 


least of quadratic order. After this assumption, the two stable coordinates can 


be written as follows: 


23 = hyyze + hy2z120 + hoz (os) 


2 2 
Z4 = $412] + $1292122 + 82225 (54) 


The coefficients h,;, s;;, of equations (53), (54), need to be determined. By 


differentiating equations (53), (54), we get the following: 


z3 = 2hy2z121 + hyo(z1Z9 =p 21 z2) + 2hooz0Zz (55) 


Za = 28412121 + $12(2122 + 2122) + 28992929 (56) 
From equation (52) and if we ignore the higher order terms, we take: 


—wz2 (57) 


24 
z2 = W004) (58) 
If we substitute equations (57), (58), into equations (55), (56), we get the 


ol 


following: 


23 = hyqwoz? a 2(ho2 — hy1)woz1z2 = hyqwoz3 (59) 
4 = 8yqwozy + 2(822 — $11)woz122 — $12wW0z5 (60) 


The third and the fourth equations of matrix equation (52) can be written in 


the following way: 
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P\z3 + [TG (Tz)](3,3) -- T~*ky3,3) (61) 


Pozq+ ee, tz eau = lime ra (62) 


z4 
In equations (61), (62), we kept terms up to second order. 


The transformation matrix T and it’s inverse T~! are 4x4 matrices, and 


their elements can be denoted by: 
i [misl, as = [roe (oad) = iL 2 3,4 (63) 


From equation (52) we have: 


dy 

TIGHT) = | 

d3 

d4 

where, 

dy = nyiG(Tz) oF moGs (Tz) 45 eke Ue Tz) (64) 
do = no G?)(Tz) + ones (Tz) als ae ‘Ge ee) (65) 
d3 = nave (Tz) ae nee (Tz) aE es: T iz) (66) 
a — niGe (Tz) + oes (Tz) + nagG? T ez) (67) 
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Also from coordinate transformation, the relationship between the old and 


new coordinates is as follows: 


p 


Q 


M4121 + M4229 + M1323 + My14Z4 


M2121 + M2222 + M323 + M7424 


M3121 + M3222 + 3323 + 3424 


M4121 + M4222 + M4323 + M4424 


(68) 
(69) 
(70) 


(71) 


Finally if we substitute equations (68), (69), (70), (71), and the expressions for 


G1, Go, G3, G4 into equations (64), (65), (66), (67), we get the final expressions 


for the coefficients d;: 


dy 


do 


dz 


d4 


= mir(liszy + laez1z2 + li7z5) 
= n19(lo5z? + log 2122 + lo729) 
9 2 
=e n13(1352} + i3gZ122 + 13729) 
= ng] (11527 aE (162129 =F l17z3) 
ale 129 (lo5 2% + 1962122 + lo7z5) 
a5 793(I35z7 a t3620220 te l37z5) 
= n3i(lyszj + 162122 + 117Z9) 
ale 739( los z4 + log2122 + lo7z5) 
+ ng33(I35z} + l362122 + 3725) 
= nai(lyszj + yezize + li7z9) 
+  m4o(losz} + logz122 + lo7z3) 


+ n43(Ig5z7 == [362122 = la7z3) 
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(73) 


(74) 


(75) 


Coefficients lj;,i = 1,2,3 7 = 5,6,7 in equations (72), (73), (74 


follows: 


Q11 2 2 
is = pH valms, + m1) 


a12 
ml 
a 


2 
TzyM3, + Iyzm3ima, + ygm21mM}1) 
L2,7731M) + ipeties + MyGM\1M3)1 


= yGW m4, ie ypBm4y) 


Q11 
li6 omg FH ¥o(2msimse + 2m, + m22) 


Q12 
+- FH (efeymsims ae Jl (m31m22 a m32™21 ) 


+ yo(maim)2 + m22m}1)| 


Q13 
FH Usy(msimee + m32M91) + 2Iy,m21m29 


+ myc(miim32 + mi2M31) + 2(ygW — yeB)maymay] 


Q11 
iz = H Yalm32 + m9) 


a12 
—(Inym3o + Iyzm3gmo2 + ygm22m 12) 


213 


COW aves ac 


(76) 


(77) 


— —[Izym32m22 + Iyzm39 + mygmigm32 + (yeW — yeB)m ia] (78) 


D 


Ne yo(m, a my) 
a2 
a D —(Ipym3, + Iyzm31mMq1 + ygmam 1) 
Q93 


2 
oe Ty aymsimar as LyzM51 + MYGM11™M3}) 


a yGW m2, = ypBm?2,) 
Q21 


log = FH yel2msimse + 2m + mo) 


a29 
5 D — |2Izy™M31M32 + Lyz(m31M2 + m32M21) 


+ ye(maimi2 + mm 1))| 
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(79) 


23 
D 
+ mycg(mi11M32 + mM12M31) + 2(yeW — yBB)maymay] (80) 


[Tay (31mg + m32M21) + 21yz2mM21M2 


G21 9 
loy = 5 ya(m3e + m2) 


A992 2 
Ae FH Ueymse + Tyz™32™22 42 YGM22™M}2) 
223 
D 
a3] 9 2 
l35 = pH yal, + M51) 


232 
me — (Iaym3y a= JUpeiekiiiton == yGM1™}1) 


[Tny™M32Ma2 + Teates + mygm2m32 + (yeW — ypB)m‘4,] (81) 


33 
D 
+ yoWmi, — ypBmj;) (82) 


; 2 
(Iry™M31M21 ae TjzM51 + mMYyGM11™M31 


231 
l3 6 — FH yel2msimse + 2m91 + m22) 


232 
ok Fp PRlaymsimse te Tyz(m31mM22 afr m32™M21 ) 


+ ye(mo21m42 + m22™M 1)| 
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_ TH [ey(maimes + m32M1) + 2lyzmMo1mM29 

+ myg(mi1m32 + m12M31) + 2(yeW — yBB)maimap] (83) 
37 = —Sya(m3, + mby) 

a —S (Leys + Iyzm32M22 + yemo22™12) 


33 
= FH Lxyims2ma2 + aoe + mygm2m32 + (yeW — ypB)m‘p] (84) 
Then equations (61), (62), can be written as follows: 
zg = 1123 + d3 t+ e3 (85) 


zg = Pozqatdgtes (86) 
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where e3, e4 come from the constant terms as: 


Ej 

Ts | 

€3 

C4 

and 

ey = kyny + konye + k3n33 
eg = kynq1 + kenge + k3n93 
e3 = kyn3, + kongq + k3ng33 
é4 = kina, + konge + kgna43 


Using equations (53), (54), we can write equations (85), (86), as follows: 
Zg = Py(hayzd + hyoz ize + hooz4) + ds + e3 (87) 
a Po(s142z7 + $1921Z9 + $9923) +d4t+e4 (88) 
And using equations (72), (73), (74), (75), they become: 
23 = (Pyhy t+ nailis + ngeles + naglss) 2} 
+ (Pyhi2 + nailig + ngaleg + 133136) z122 
+ (Pyhog + nailiz + nggle7 + nggls7)z5 + €3 (89) 
ég = (Posi. + nailis + ngglos + n33l35) 24 
+ (P2812 + nagilig + ngalag + 233136) 2122 
+ (P2829 + nailiz + ngele7 + naglg7)z5 + e4 (90) 
Comparing the coefficients of equations (89), (90), with the coefficients of 
equations (59), (60), we get: 
—Pyhy, + wohi2 = 131115 + gales + 133135 
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—2wohi — Pihi2 + 2woha2 = n3ilig + N3aleg + n33l36 


—wohi2 — Pyhog = 231117 + N32l07 + 233137 


Solution of the above 3x3 linear system of equations, gives us coefficients hj, 


hy2, hoo. 
Also in the same way: 
— P9811 + w0si2 = Naylis + N4gle5 + 243135 


—2wo9811 — P2812 + 2wos22 = Nailig + Nagleg + 43136 


—W9812 — Po822 = naili7 + nagle7 + 243137 


Solution of the above 3x3 linear system of equations, gives us coefficients 511, 


$12, 522- 


E. AVERAGING 


In this part of our analysis we are going to take into account the third order 


terms of equation (52): 


ii 
Pei A 
I 
Where, 
fi = my G? (Pz) + nieGo(3)\(Tz) + nigGh? (Tz) (91) 
fo = nynG? (Tz) + noeGo(3)(Tz) + n3xG$? (Tz) (92) 
fg = ngiG2 (Pz) + naGo(3\(Tz) +nasGs (Tz) (93) 
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3 3 
fa = na GP(Tz) + nseGo(3)(Tz) + ngGP(Tz) (94) 
If we substitute equations (68), (69), (70), (71) and the expressions for Gj, 
G2, G3 into equations (91) throuh (94), we get: 
_ 3 2 D 3 
fi = 34 (l41 2} ar l192} Vay) =P 1132125 == 11425) 
+ n42(lo1z9 + log z? z9 > logzizs oF loazs) 
3 2 2 3 
A n43(l312} + l39z3 22 + 1332125 aa I34z5 ) (95) 
fo = moy(laiz? + lezeze + lygzizs + ly4z3) 
ar n99(lo1z° == logz? 29 == lo3z125 + loz) 
3 2 2 3 
+ 93 (1312} + 13923 22 = 1332125 + 13425) (96) 
f3 = mgi(liiz; =e lygz? 29 =e 432125 Sp lyaz>) 
a6 n39(lo12° + logz? zo + logzize -- logz2) 
3 2 2 3 
te 733 (13124 + 1392729 + 1332129 + 13425 ) (97) 
f4 = nay(14422 + lygz? zo + lygz123 + L4z5) 
ar n4o(lo12z3 Se log z zo a logz1z5 AF logzs) 
3 2 2 3 
ae n43(131z} oF l30Z4 Vig) a 1332125 ae 13425) (98) 
From equation (52) and from the system matrix, the derivatives of the first 


two modified coordinates, z; and Zz. can be written as follows: 


Zz = (a’€)z] — (wo twe)zo + FF (21, 22) (99) 
zg = (wotw'e)zy +@ezq + FFo(z1, 2) (100) 

where, 
Pg Bo Sn sen ae ea (101) 
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FFo(z1,22) = dot fetes (102) 


If we combine equations (101) and (102) with equations (72) through (75), 


and (95) through (98), we get: 


FF (21, 29) — r42z> oa BUSES a 1132124 ae 11425 
+ pyu2? t+ piezize + pi3z5 + €1 (103) 
a. a Z 2 5 

F F9(21, 29) = 1912) + 1222122 + 1232125 + 1425 
+ poizy + pooz1z2 + pogz5 + 1 (104) 


where coefficients r;; and p,; are: 


Ti = Nl + Ni2Ql21 + 213/31 
Ty2 = Nil + NygQleq + N1i3l32 
713 = 111113 + Niele3 + 213133 
Tia = Nyl14 + Nyele4 + Niglza 
rai = nailiy + naglei + no3l31 
Too = Nailig + Noagloe + N3139 
T93 = N21113 + Noaglo3 + n3133 


Toa = Noili4 + Noglog + Nagla 
or generally, 
Tig = Mah; + Nolo; + niglz, t= 1,2 7 =1,2,3,4 (105) 
also, 


Pir = NMiylis + Ny2le5 + 213/35 
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Piz = Nile + N12lo6 + n13ls¢ 
P13 = M117 + Ny2le7 + 133137 
pai = Nails + nogles + no3lgs 
P22 = Nailig + Noalog + no3ls¢ 


P23 = N21117 + Noegle7 + n93137 


or generally, 


Pig = Malik + Niglee + gl, t= 1,2 g=1,2,3 kK=j7+4 (106) 


The coefficients 1,; 1=1,2,3 j =1,2,3,4 are as follows: 


liy 


lio 


_ au Cd, 

FOB os 
1 
Ale = B)mi1] 
ai2 Cp, 


D' 6y 


1 a 
5 (teW ~ 2pB)m3,] + ap (eeW — zpB)m3, (107) 





(Egm?, a 3Eym2,m21 ae 3E9m ym, oe E3m3,) 





(Eym?, + 3Eom?,ma + 3E3my1m3, + E4m3,) 


411.-Cp 
— BGP BBominmis == df (m2,ma9 -- 2m11™42M21) 
1 
3Eo(myms, =e 2m1M22M}1) a 3E3m3,m22) ae AC ~ B)3m2,ma9] 
012 Dy | 
DD 67 


3E3((mio9m3, == 2m21M22™M}1 ) a 3E4m3,mo) 


3Eymi m1 + 3Eo((m2,m2 + 2m44m 2™M1) 


> (zaW — z—B)3m 3, map (108) 


(3E9mym?, + 31 (m2ymay a 2m 11M 42M22) 


1 
5 taW = tp B)3m2,ma9] lis 


_ a1 Edy 
DY 67 





3Eo(m5omiy + 2mo1m22M12) + 3E3mo1M5p) 
1 


rAd — B)3maim4o} 
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lia 


loi 


loo 


lo3 


a12 ,CDy 
D 6-7 


3E3(moomi1 6 2m21M22™M12) a 3E4mo1m5o) 


——+(3Eym11m?, + 3Eo(m2,m2 + 2M 141M12™M2) 


i Q13 
=(tGW — 2—B)3mamip] + a (zeW — zpB)3mamig 





211 Cd, 
D~ 6y 


=(W a B)m4p] 


(Eym3, + 3E\m? may + 3E2m12Ms5 + E3m3,) 


a12 CD, 
D~ 67 


1 
5 \zaW — tpB)m3,] + 


(Eym3, ae 3Eom<)ma ae 3E3m42mMs5 94 E4m3,) 





7 el — zpB)m3, 


ao Cp 
——=[—* Eom’, + 3Eymi,mo1 + 3E2m 41m, + E3m3,) 


Gay, 
1 ‘ 
iW _ Bm, 





22 CD, 
D~ 6y 


it 
g(zaW = tpB)\m3,|+ 


(Eym?, + 3Eqm?,m2 -+- 3E3m41m3, + E4m3,) 





= —(z0W = zpB)mi, 


mee C'p 
: = |—— 2 (3Egmi,mi2 + 36} (m2 mo + 2mymM12M1) 


D~* 6y 


(109) 


(110) 


(111) 


1 
3Eo(mypms, I= 2M21M22M}41) =e 3E3m31m2) + 5 iW — B)3m3,m49] 


22 (CDy 
Deon 


3E3((myoms, ss 2™21M22M}41) ate 3E4m3s,m22) 





(3Eym?,mo + 3Eo((m?,mo2 + 2m41m12M1) 


I 
gizaW — tpB)3mi,mao] + ——(zgW — zpB)3mi,m42 


5 


a C 
a —|(—+ Py (3E9m1m4p + 3Fj (m2ym + 2M 11™M12M22) 


D~* 6y 
3Eo(msomy1 + 2mgymM12) + 3E3mo1M5q) 


] 
A =a B)3m4,m4p] 


a22 Cp 
= ey (Z3Eymymjg + 3E2(migma; + 2mymgM29) 





Al 


Qa) 
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3E3(m5gmi + 2mgimggm12) + 3E4mo1M5q) 


] 


; se (eG = zpB)3m41mig 


—(rgW = tp B)3m4,m%,] Siete 6D 


221 C'p, 
~(Eo 9m. + 3Eym4,m + 3Egmgmin + E3m3o) 


D Dey oy 
(WW = B)m io] 


ao .C'p 
Dl oy im ai 3Eom* ma + 3E3m12M55 ag E4m)) 
= (tGW a tpB)m%,] ig ae _ zpB)mio 


431 31 /CDy 


D ay (Eom?, == 3Eym~)ma == 3Egm14m3, =s E3m3,) 


=(W zr B)m4,] 





032 ae? 
De Gay 


1 
5 \toW = tpB)m3,] + 





11 + 3Egmjymai + 3E3m11m3, + Esm};) 


= —=(zcW _ zpB)m}, 


_ 431 a31/C Dy 


D by (3E9m2m1 Se df (m?,mo2 == 211M 42M) 


(113) 


(114) 


(115) 


2 1 2 
3E2(mi9m5, ala 2m21M992™M}1) + 3£3m51M22) aR rAd = B)3m4,mMa9] 


a3q.C 
3E3((myom3, == alae + 3E4m31m29) 


= (GW — rpB)3m‘i,mao] is ap (zaW — zpB)3mi,;ma2 


"ani Co, 

DS 67 
3Eo(msomi1 + 2mo1™m22™12) + 3E3m21mM5q) 
1 


A = B)3ma4ymio] 


932 (© Dy 
D6 


3E3(meomiy + 2mo1m22™m12) + 3E4m21m5p) 


(3Eom11m?2, + 3F} — as 241M 12™M22) 


(3E;mym?, == 3E2(m?.mo ao 2m41™m12M29) 
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(116) 


a 


1 a 
5 \teW = t BB)3m4\m‘4o| AR ap 2aW _ 2pB)3maim4p 


a31,Cp 
—| - (Eom, = 3Eym?,ma1 =F 3E gm ym, 5= E3m35) 


 D* by 
1 
5 \W im B) mo] 








032 /CDy 
D~ 6y 


if 233 
g zaW — x pB)mo| + ap ew — zpB)m3, 


(E,m3, a 3Eomi4ma + 3E3m12M55 + E4m3p) 


The next step is to introduce polar coordinates in the form: 


24, > Rceos@ 


22> Resin @ 


(117) 


(118) 


(119) 


(120) 


We use polar coordinates, because it is easier in this way to investigate the 


existence of limit cycles. 


Substituting equations (119), (120), into equations (99), (100), we get: 


Rcos@ — R(sin@)@ = (wo+w'e)Rcosé +a'eRsind + 


P,(0)R° + Qi(0)R? +e, (121) 
Rsin@+R(cos0)@ = (wot+w'e)Rcosd +a'eRsind + 
P2(8)R° + Q2(0)R* + €2 (122) 
where, 

Ce — ea cos? @ + rip cos” 6 sin 6 + 113 cos 6 sin? 6+ 7ri4 sin? @ (128)) 
On) on cos? 6 + pip cos 4 sin @ + p33 sin? 6 (124) 
EGC) — 55 cos® @ + reg cos” @ sin 8 + 793 cos 4 sin? 6 + roq sin°@ (125) 
COD) ta cos” 6 + po cos @ sin 6 + pog sin? 6 (126) 
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If we multiply equation (121) by cos@ and equation (122) by sin@, and add 


the two resulting equations, we get: 


R = a'eR+ P(0)R* + Q(0)R? + (e; cos 6 + eg sin 6) (127) 

where, 
P(6) = P,(@)cos@ + P9(6@) siné (128) 
Q(@) = Q1(0) cos 6 + Q2(0) sin 8 (129) 


Equation (127) contains one variable that varies slowly in time (R) and a fast 


variable (@). 


If we average this equation over one complete cycle in 0, from 0 to 27, 


equation (127) becomes: 


R=a'eR+KR°3+LR?+M (130) 
where, 

1 PA 9 9 

a — | P(6)d 
27 JO ( ) 
1 

= g(oru + 113 + 722 + 3124 (131) 

1 Pe 

ia — | Q(6) dé (132) 
27 JO 


ih 27 
1 — | e€, cos 6 dé 
27 JO 
1 2r 
+-_ — / e5 sin 6 dé 
27 JO 


= 5~[sin — =~ [cos ae 


0 (133) 
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Finally equation (135) becomes: 


R=a'eR+KR?® (134) 


F. LIMIT CYCLE ANALYSIS 


At steady state R = 0,and equation (134) becomes: 
0 = R(a'e + KR’) (135) 


Equation (135) has two solutions. The first solution is R = 0. This is the 
trivial solution and it does not give us much information.The second solution 


1S: 
—ale 


Te 





pe (136) 


This solution gives us a limit cycle of constant amplitude R in the 21, z 
cartesian coordinate system.This limit cycle exists if the quantity inside the 


square root is possitive, or 
—a’e 


K 





> 0 (137) 


Condition (137) is necessary for the amplitude of the limit cycle, R, to be a 


real number. 


In our case a’ is always negative, because for constant rg, as we decrease 
« (Figure 1), the real part of the critical pair of eigenvalues increases, due to 


further loss of stability. In other words we can say that: 


a a0 (138) 
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From conditions (137), (138) we see that the existence of the limit cycles 


depends on the value of parameter kK. We can see that: 


e Ii kK < 0; periodic solutions exist for « < 0 or zg — zg, < 0 or zg < Zzg., 


and they are stable. 


e Ii K > QO, periodic solutions exist for « > 0 or zg ~ zg, > 0 or zg > 2g,, 


and they are unstable. 


The characteristic root of equation (134) in the vicinity of (136) is: 
B= —2a’e (139) 


The sign of this characteristic root, assigns the stability of the periodic solu- 


tions. 


G. RESULTS AND DISCUSSION 


Typical results in terms of the nonlinear stability coefficient K are pre- 
sented in Figures 3 through 6. The stability coefficient K is shown in its 
normalized form as K - y (Papadimitriou, 1994). Figure 3 shows K versus 
the LCG/LCB separation distance xg (in ft) for a given vehicle speed and for 
different values of the drag coefficient. It can be seen that K is everywhere 
negative, which means that all bifurcations to periodic solutions are super- 
critical. Higher values of the drag coefficient result in stronger supercritical 
bifurcations which means that the corresponding limit cycle amplitudes will 


be smaller. 
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Figure 3: K - y versus xg for U = 5 ft/sec and different values of Cp, 





-6 


XG 


Figure 4: K -y versus xg for Cp, = 0.5 and different values of U (ft/sec) 
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Figure 5: Simulation results (¢,t) for Cp, = 0.5, U = 5 ft/sec, and zg = 1 ft 
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Figure 6: Limit cycle amplitudes (in ¢) versus zg for U = 5 ft/sec and rg = 1 ft 
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. 
a 
a 





Figure 4 shows K versus zg for a given drag coefficient and for different 
forward speeds. It can be seen that the bifurcations are supercritical with the 
possible exception of large speeds and small zg. This means that it is possible 
for a properly trimmed vehicle at relatively high speeds to experience an oscil- 
latory behavior even before stability is lost. This demonstrates a destabilizing 


effect which could not have been predicted by linear techniques. 


Figures 5 and 6 shows numerical simulation results which confirm the theo- 
retical predictions. Both figures correspond to a supercritical bifurcation case 
and they show a continuous increase in limit cycle amplitudes as the bifurca- 


tion point is crossed. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


This thesis presented a comprehensive nonlinear study of straight line sta- 
bility of motion of submersibles in coupled sway/yaw/roll motions under open 
loop conditions. Primary loss of stability was shown to occur in the form of 
Hopf bifurcations to periodic solutions. This loss of stability is characteristic 
of the coupling of roll into sway and yaw and cannot be predicted by consid- 
ering the uncoupled motions. The critical point where instability occurs was 
computed in terms of vehicle metacentric height, longitudinal separation of 
the centers of buoyancy and gravity, and the forward speed. Analysis of the 
periodic solutions that resulted from the Hopf bifurcations was accomplished 
through Taylor expansions, up to third order, of the equations of motion. A 
consistent approximation, utilizing the generalized gradient, was used to study 
the non~—analytic quadratic cross flow integral drag terms. The results indi- 
cated that loss of stability occurs always in the form of supercritical Hopf 
bifurcations with stable limit cycles. It was shown that this is mainly due to 
the stabilizing effect of the drag forces at high angles of attack. Subcritical 
bifurcations, however, with considerably higher limit cycle amplitudes may 
develop for sufficiently high forward speeds and small LCG/LCB separations. 
Simulation studies of these subcritical bifurcations along with the effects of 


vertical plane coupling constitute recommendations for further research. 


ol 





APPENDIX 


The following is a list and description of the computer programs used in this 
thesis. The programs are written in FORTRAN or MATLAB. Complete print- 


outs of the programs follow after the list. 


e STABILITY.M 


MATLAB program for performing linear stability analysis. 


e SIM.M 


MATLAB program and functions for numerical simulation. 


e HOPF.FOR 
FORTRAN program for calculating the nonlinear stability analysis coef- 
ficient K. It requires data from STABILITY.M and standard numerical 


linear algebra subroutines. 


o3 


FeorAb ele ry el 

h 

7% LOSS OF STABILITY 

i, HK KK eK KK KK OK KK KK KOK OK KK KK OK KK KK OK KK KK KK KKK KKK 
a=1 

W=12000; 

IXX=1760; IYY=9450; 

IZZ=10700; IXZ=0; IXY=0; 

TYZ=0; L=17.425; RHO=1.94; 

G=32.2; U=6.5; M=W/G; B=W; 


ND1=0 .5*RHO*L™2; 
/% DEFINE HYDRODYNAMIC COEFFICIENTS 


YPDOT=1 .270e-04*ND1*L“2; 
YVDOT=-5 .550e-02*ND1*L; 
YRDOT=1 . 240e-03*ND1*L°2; 
YP=3 .055e-03*ND1*L; 
YV=-9.310e-02*ND1; 

YR=-5 .940e-02*ND1*L; 


NPDOT=—-3.370e-05*ND1*L°3; 
NVDOT=1 .240e-03*ND1*L~ 2; 
NRDOT=-3. 400e-03*ND1*L73; 
NP=-8 .405e-04*ND1*L" 2; 
NV=-1.484e-02*ND1*L; 
NR=-1 .640e-02*ND1*L*2; 


KPDOT=-1 .01e-03*ND1*L73; 
KVDOT=1 .27e-04*ND1*L“ 2; 
KRDOT=-3.3/7e-05*ND1*L73; 
KP=-1 .10e-02*ND1*L*~2; 
KV=3 .055e-03*ND1*L; 
KR=-8 .41e-04*ND1*L*2; 


flag=0; 
for XG=0:0.01:2, 


flag=flag+1; 


o4 


xg (flag)=XG; 

a=IXX-KPDOT; b=KP*U; e=KV*U; 
f=KRDOT; i=YP*U; j=M-YVDOT; k=YV+*U; 
1=XG*M-YRDOT; m=U*(CYR-M); o=NPDOT; 
p=NP+U; q=-XG*W; r=XG*M-NVDOT; 

w=U* (NR-XG*M) ; x=NV*U; u=IZZ-NRDOT; 


al=-u*M"2; 
a2=—-u*M* YPDOT—-u*M*KVDOT+M*0*1+f£%*r*M ; 
a3=a*j*xu-a*l*r—u*KVDOT*YPDOT+KVDOT*0*1+£*r*YPDOT-f£*0*]j ; 


b1= (wx (M72) )+(r*U* (M72) ) ; 
b2=-M*e*u-M*u*i+w*M* YPDOT+w*KVDOT*M—o*m*M+p*1*M-£*x*M+. .. 
r*M*U*YPDOT-—o* j *M*#U+r*KR*U*M ; 


b3=-e*u* YPDOT+e*o*1-u*i*KVDOT+w*KVDOT*YPDOT-KVDOT*o*m+KVDOT*p*1-... 


ax j*w-a*xk*euta*] *xt+a*xrem—b*j*utb*l*r+f*r*i-f£*x*YPDOT+... 
o*k*f£-f*p* j+r*KR*U*YPDOT ; 


c1=-x* (M72) *U; 
C2=r*i*M*U-x*M*U*YPDOT-x*KR*U*M+0*k*M*U—p* j *M*U+j*u*xW-lL*r+*w-... 
qx 1l*Mt+w*i*M-m*p*M+e*wn ; 
c3=a*xk*w-axm*xtb* ] *w+b*k*u-b*1l*¥x—-b¥r*m+r*i*KR*U-x*KR*U*YPDOT+... 


o*k*KR*U-px* j *KR*¥U-g*1*KVDOT+w*i*KVDOT-m*p*KVDOT-e*u*ite*w*YPDOT-... 


exo*mte*p*]l+f£*q* j-L£*x*i+f*pek 5 


d2=q* ]*M*U-x¥*i*M*U+p*k*M*U+q4m*M—j *eweW-k*utwt+]l*x*Wt+rem 5 
d3=q*j*KR*U-x*i*KR*U+p*k*KR*U-f*q*k+qem*KVDOT—-exg*l+e*wei-... 
e*m*p—b*k*w+b*m*x ; 


e2=k*w*W-m*x*W-g*k*Mx ; 
e3=e*q*m—q*ek*KR*U 5 


£5=(-e2* (b172))+b1i*c1*d2; 

£4=( (b2*c1+b1*c2) *d2)+b1*c1*d3-al* (d2°2) -e3* (b17 2) -2*e2*b1*b?2; 
f3=-a2* (d272)-2*d3*d2*a1-2*e3*b1*b2—-—e2%* ((b272)+2*b1*b3) +... 
d2* (b3*c1+b2*c2+b1*c3)+d3* (b2*c1+b1¥*c2) ; 

f2=-e3* ((b2°72)+2*b1*b3) -2*e2*b2*b3+d2* (b3*C2+b2*C3) +... 

d3* (b3*c1+b2*c2+b1*c3)—-a3* (d2~2)-2*d3*d2*a2—-al* (d3°2) ; 

f 1=b3*c3*d2+d3* (b3*C2+b2*c3) -2*e3*b2*b3-e2* (b372)-... 
2*d3*d2*a3-a2* (d3~2) ; 

£0=b3*c3*d3-a3* (d3~2) -e3* (b372) ; 


oO 


eoct= oem t5 922 £1 £0)’. 
ZG=uaets (coef) : 


tot (flag)=ZG(5,1); 
end 


plomexe, Lot)4prid; 

title(’ZG versus XG plot in the point of loss of stability’) 
xlabel (7 xG. inet’); 

ylabel(’ZG in ft’); 


06 


% NON LINEAR SIMULATION PROGRAM 


tO=0; 

tfinal=500; 

q=(1/180) *pi; 

yo=[0 0 0 q]; 

[t, y]=ode45(’vdpo2’ ,t0,tfinal,y0) ; 


fupure(1) ,.plot(t(:,1),57.29578*y(: ,4)) ,grid 


title(’plot of fi with time’); 
ylabel(’fi’) ,xlabel(’time’ ) 


of 


% NON LINEAR SIMULATION PROGRAM 


function yprime=vdpo2(t,y) 
W=12000; 

IXX=1760; IYY=9450; 

IZZ=10700; IXZ=0; IXY=0; 

TYZ=0; L=17.425; RHO=1.94; ZB=0; 
G=32.2; U=5; M=W/G; B=W; XB=0; 
ZG=0.05; CD=0.5; YG=0; 

XG=1; YB=0; 


ND1=0.5*RHO*L™ 2; 
% DEFINE HYDRODYNAMIC COEFFICIENTS 


YPDOT=1.270e-04*ND1*L* 2; 
YVDOT=-5 .550e-02*ND1*L; 
YRDOT=1 .240e-03*ND1*L*2; 
YP=3.055e-03*ND1#L; 
YV=-9.310e-O2*ND1 ; 

YR=-5 .940e-O02*ND1*L; 


NPDOT=-3.370e-05*ND1*L73; 
NVDOT=1 .240e-03*ND1*L* 2; 
NRDOT=-3.400e-03*ND1*L73; 
NP=-8 .405e-O04*ND1*L°2; 
NV=-1.484e-02*ND1*L; 
NR=-1.640e-02*ND1*L° 2; 


KPDOT=-1 .O1e-03*ND1*L73; 
KVDOT=1 .27e-O04*ND1*L~2; 
KRDOT=-3.37e-05*ND1*L* 3; 
KP=St LOC-OZ*NDI+L 2; 
KV=3.055e-03*ND1#L ; 
KR=-8.41e-O4*ND1i*L° 2; 


D=( (M-YVDOT) * (IZZ-NRDOT) * (IXX-KPDOT) )... 
—((M-YVDOT) * (IXZ-KRDOT) * (-IXZ-NPDOT))... 

— ((M*XG-NVDOT) * (M*XG-YRDOT) * (IXX-KPDOT) )... 

+ ((M*XG-NVDOT) * (IXZ-KRDOT) * (-M*ZG-YPDOT))... 
+ ((-M*ZG-KVDOT) * (M*XG-YRDOT) * (-IXZ-NPDOT))... 


—((-M*ZG-KVDOT) * (IZZ-NRDOT) « (~M*ZG-YPDOT) ) ; 


A11=( (IZZ-NRDOT) * (IXX-KPDOT) ) - ((IXZ-KRDOT) * (-IXZ-NPDOT)) ; 
A12=( (-M*XG+YRDOT) * (IXX-KPDOT) )+( (IXZ-KRDOT) * (-M*ZG-YPDOT)) : 
A13=( (M*XG-YRDOT) * (-IXZ-NPDOT) )-( (IZZ-NRDOT) * (~-M*ZG-YPDOT)) ; 
A21=((-M*XG+NVDOT) * (IXX-KPDOT) )+( (-M*ZG-KVDOT) * (-IXZ-NPDOT)) ; 
A22=((M-YVDOT) * (IXX-KPDOT) ) -( (-M*ZG-KVDOT) * (-M*ZG-YPDOT)) ; 
A23=((-M+YVDOT) * (-IXZ-NPDOT) ) +( (M*XG-NVDOT) * (-M*ZG-YPDOT) ) ; 
A31=( (M*XG-NVDOT) * (IXZ-KRDOT) ) ~ ((-M*ZG-KVDOT) * (IZZ-NRDOT)) ; 
A32=( (-M+YVDOT) * (IXZ-KRDOT) ) + ( (-M*ZG-KVDOT) * (M*XG-YRDOT) ) ; 
A33=((M-YVDOT) * (IZZ-NRDOT) )- ( (M*XG-NVDOT) * (M*XG-YRDOT)) ; 


ie EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 

h 

F(1,1)=(A11*#YV*U+A12*NV*U+A13*KV4U)/D; 

F(1,2)=(A11* (YR*U-M*U) +A12* (—-M*XG*U+NR*U) +A13* (M*ZG*U+KRXU) ) /D; 
F(1,3)=(A11*YP*U+A124*NP*U+A13*KP*U) /D; 

F(1,4)=0; 

F(2,1)=(A21*YV*U+A22*NV*U+A23*KV4*U) /D; 

F(2,2)=(A21* (YR*U-M*U) +A22* (-M*XG*U+NR4*U) +A23% (M*ZG*U+KR*U) ) /D; 
F(2,3)=(A21*YP*U+A22*NP*U+A23*KP*U) /D; 

F(2,4)=0; 

F (3, 1)=(A31*YV*U+A324NV*U+A33*KV4U) /D; 

F(3,2)=(A31* (YR*U-M4#U) +A32* (—M*XG*U+NR*U) +A33* (M*ZG*U+KRU) ) /D; 
F(3,3)=(A31*YP*U+A32*NP*U+A33*KP*U) /D; 

F(3,4)=0; 

F(4,1)=0; 

F(4,2)=0; 

F(4,3)=1; 


% DEFINE THE LENGTH X AND HEIGHT H TERMS FOR THE INTEGRATION 
poe INerEET . 

REC1)=-105- 97 12; 

XL(2)=-104.3/12; 

RECS) ==99).3/7 12; 

XL(4)=-94.3/12; 

HE (S)=-87.3/ 12: 

XL(6)=-76.8/12; 
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XL(7)=-66.3/12; 
XL(8)=-55.8/12; 
KL(9)=72.7/ W2s 

ML (10) =79.42/7 12: 
XLCi1) =6322/ 12. 
XL(12)=87 .2/12; 
ALCS)H=91 227 i 
XL(14)=95.2/12; 
XL(15)=99.2/12; 
KEC16)=101-27 12; 
XLC17)=102. 1712: 
RECS =103-2/ 12: 


HT (1) =0; 

Hr (2223712. 
HT (3) =8.24/12; 
HT (4) =13.96/12; 
HT (5)=19.76/12; 
HT(6)=25.1/12; 
HT (7) =29.36/12; 
HT (8) =31.85/12; 
HT (9) =31.85/12; 
HT (10)=30.00/12; 
HT (11)=27 .84/12; 
HP(12)=25.12/12; 
HT (13) =21.44/12; 
HT (14)=17.12/12; 
HT (15) =12.0/12; 
HT(16)=9.12/12; 
HT (17) =6.72/12; 
HT (18)=0; 


bone kal wiltsi 

VEC1 (i) =HT Ci) * Cy (1) +XL (i) *y (2) ) * Cabs Cy (1) +XL Ci) *y (2) )) ; 
VEC2(i)=HT Ci) * Cy (1) +XL (i) *y (2) ) * Cabs (Cy (1) #XL (i) *y (2) )) *XL(i) ; 
end 


OUT=0; 

for jE: i, 

OQUTI=0.5¥ (VEGI (7) )+VEGI (+1) (XL(j) +1) -XL(j)); 
OUT=OUT+0UT1; 

end 
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IV=CD*OQUT ; 


OUT=0; 

Lor j= Lald, 

OUT2=0 .5* (VEC2(j )+VEC2(j+1) )* (XL(jt+1) -XL(j)); 
OUT=OUT+QUT2; 

end 

IR=CD*QUT ; 


yprime=[F(1,1)*y(1)+F(1,2)*y(2)+F(1,3)*y(3)+F(1,4)*y(4)+... 
(A11/D) * (YG* Cy (3) 72+y (2) ~2) -IV+(W-B) *sin(y(4)))+... 


(A12/D) * (IXY*y (3) “2+1YZ*y (3) *y (2) +YG*y (1) *y (2) -IR+ (XG#W-XB*B) *sin(y(4)))+... 


(A13/D) * (-IXY*y (3) *y (2) -LYZ*y (2) “2-M*YG*y (1) *y (3) +(YG*W-YB*B) *cos(y(4))... 
-(ZG*W-ZB*B) *sin(y(4)));... 

F(2,1)*y(1)+F(2,2) *y(2)+F(2,3) *y(3)+F (2,4) *y(4)+... 
(A21/D) * (YG* (Cy (3) “2+y (2) 72) -IV+(W-B) *sin(y(4)))+... 
(A22/D) * (IXY*y (3) “2+1YZ*y (3) *y (2) +YG*y (1) #y (2) -IR+ (XG*W-XB*B) ¥sin(y(4)))+.. 
(A23/D) * (-IXY*y (3) #y (2) -LYZ¥y (2) *2-M#*Y Gy (1) ty (3) + (YG*W-YB*B) *cos(y(4))... 
- (ZG*W-ZB*B) *sin(y(4)));... 

F(3,1)*y(1)+F (3,2) ¥*y (2) +F (3,3) *y(3)+F (3,4) *y (4)+... 
(A31/D) * (YG* Cy (3) *2+y(2)~2)-I1V+(W-B) *sin(y(4)))+... 


(A32/D) * (IXY*y (3) 72+1YZ*y (3) *y (2) +YG*y (1) *y (2) -IR+ (XG#W-XB*B) *sin(y(4)))+... 


(A33/D) *(-IXY*y (3) *y (2) -LYZ*y (2) ~2-M*YG+y (1) *y (3) + (YG*W-YB*B) *cos(y(4))... 
—(ZG*W-ZB*B)*sin(y(4)));... 
1*y(3)]; 
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Ci) Gane) 


HOPE POR 


HOPE 5 


IMPELIC 
REAL*8 
REAL*8 
REAL*8 
REAL*8 


REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 


IFURCATIONS PROGRAM 


ii DOUBLE PRECIS@ON CA-H,0-Z) 

Peny oie YPDOT, YVDOT {NDI YRDOT 
YP,YV,YR,NPDOT,NVDOT,NRDOT,NP,NV,NR 
GAMA,U,KVDOT,KRDOT ,KPDOT,D 

BO El E2pes £45 3G,Z2Gynk, KPPKVeNG! ,ZG1 


Pe ieee la i2i M22 M23 

M24 ,M31,M32,M33,M34,M41,M42,M43,M44 
N1i1,N12,N13,N14,N21,N22 ,N23 ,N24 
N31,N32,N33,N34,N41,N42,N43,N44 

Bit vi2 bie skl4 L221 ,L22,L23,L24,1L31 

L32 ,L33 ,L34,L15,L16,L17 ,L25,L26,L27,L35 
L36,L37,L1A,L2A,L3A,L4A,L5A,L6A,L7A,L8A,L9A 
DIORA LIGA, LIi2a L122, L3,L4,05,L6,L7 

Be 709 RIT R12, R13, R14,R21,R22,R23, R24 

Pam Pi2 Plo, P21,P22,P23 


DIMENSION F(4,4) ,T(4,4) , TINV(4,4) ,FV1(4) ,IV1(4) ,YYY(4, 4) 


DIMENS 
DIMENS 
DIMENS 


ION WR(4) ,WI(4) , TSAVE(4,4) , TLUD(4,4) , IVLUD(4) , SVLUD (4) 
ION ASAVE(4,4) ,A2(4,4) ,XL(18) ,HT(18) , ZGG(197) , FF(4,4) 
ION VECO(18) , VEC1 (18) , VEC2(18) , VEC3(18) , VEC4(18) ,XGG(197) 


NGG ieee.) 


OPEN ( 
OPEN ( 
OPEN ( 
OPEN ( 


WEIGHT 
IXX 
Ig 6 
TAZ 
IXZ 
TeXGY 
YZ 

Ik 

RHO 


20, FILE=’HOPF.RES’ ,STATUS=’ OLD’ ) 
21,FILE=’DATA.DAT’ ,STATUS=’ OLD’ ) 
23 ,FILE=’HOPF1.RES’ ,STATUS=’ OLD’ ) 
25, PILE="HOPF2-RES” ,STATUS=" OLD’ ) 


=12000.0 
=1760.0 
=9450.0 
=10700.0 
=0.0 
=0.0 
=0.0 
=17.425 
=1.94 
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@ 


a 


WRITE (+,4) > ENTER CD’ 
READ (*,*) CD 
G =32.2 

XB =0.0 

ZB =0.0 

YG =0.0 

YB =0.0 
YDELTAR=0.0 
DELTAR=0.0 
NDELTAR=0.0 
NPROP=0.0 

M= WEIGHT/G 

B= WEIGHT 
W=WEIGHT 

U=8.0 
ND1=0.5*RHO*L**2 


i 


DEFINE HYDRODYNAMIC COEFFICIENTS 


YPDOT=1 .270E-04*ND1*L**2 
YVDOT=—5 .. 550E-02*ND1*L 
YRDOT=1 . 240E-03*ND1 *L**2 
YP=3.O55E-03*ND1i+L 
YV=-9.310E-O2*ND1 

YR=-5 . 940E-02*ND1*L 


NPDOT=-3.370E-05*ND1*L*¥*3 
NVDOT=1 . 240E-03*ND1*L**2 
NRDOT=-3 . 400E-03*ND1*L**3 
NP=-8. 405E-04*ND1*L**2 
NV=-1.484E-02*ND1*L 

NR=-1 .640E-02*ND1i*L**2 


KPDOT=-1.01E-03*ND1*L**3 
KVDOT=1 .27E-04*ND1*L**2 
KRDOT=~-3 .37E-O5*ND1*L**3 
KP=-1 . 10OE-O2*ND1*L**2 
KV=3.O055E-03*ND1*L 

KR=-8 .41E-O04*ND1*L**2 


DEFINE THE LENGTH X AND HEIGHT H TERMS FOR 
THE INTEGRATIONS SALE IN FEET: 
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ME Gel) =—106,,9712.0 


er ¢ 
XL( 
XL 
XL( 
XL( 
EL 
XL( 
KL 


2)=-10 
3)=-99 


4)=-94. 


5) =-87 
6)=-76 
7) =-66 
8)=-55 
9)=72. 


XL(10)=79. 
XL(11)=83. 
XL(12)=87. 
XL(13)=91. 
XL(14)=95. 
XL(15)=99. 
XL(16)=101.2/12.0 
XL(17)=102.1/12.0 
KECie)—=10522/712..0 


HT( 1)= 0. 
BiG) — 2) 
HiGs)= Ss. 
HT( 4)= 13. 
HT( 5)= 19 
Hr@o)= 25. 
HT( 7)= 29 
HT( 8)= 31 
HT( 9)= 31 
HT(10)= 30 
HT(11)= 27 
HT(12)= 25. 
HT(13)= 21 
HT(14)= 17. 
HT(15)= 12 
HT(16)= 9. 
HT(17)= 6. 
HT(18)= 0. 


4.3/12.0 
07 ler 
3/12. 
53/126 
poy LP) 
£o/ U2 
nO7 12. 
G7 2. 
27 ze 
PL) SO 
2/12. 
Py uPA 
Dye 
Payee 


OO OOo © 


OOO 0 0 0 OO 


000 
28/1220 
24/12.0 
3672-0 
SiS ARO) 
1/12.0 
T6712. 
p50/ 12), 
souy 12. 
B007 125 
84/12. 
eZee 
44/12. 
2/12. 
707 12-0 
127 120 
“2712.0 
00 


O CO OFC. o 2 © 


DO 104 K = 1,18 
VECO (K) =HT (K) 


VEC1i (K)=XL (K) *HT (K) 
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Q 


Q 


104 


& 


Rk & & 


VEC2(K) =XL (K) *XL(K) *HT (K) 
VEC3(K)=XL (K) *XL (K) *XL (K) *HT (K) 
VEC4 (K) =XL (K) *XL(K) *XL(K) *XL (K) *HT (K) 
CONTINUE 
CALL TRAP(18,VECO,XL,EO) 
CALL TRAP(18,VEC1,XL,E1) 
CALL TRAP(18,VEC2,XL,E2) 
CALL TRAP(18,VEC3,XL,E3) 
CALL TRAP(18, VEC4,XL,E4) 


WRITE (*,*) ’ ENTER GAMA’ 
READ (*,*) GAMA 


ee ce ce ce cr i ce rr cr cr cc cr cc rr re rc re ee ee ee ee ee 
SS cep ce erp ces re re ec cc ee er a ce ce ce ce ee ee ee ee ee ee ee ee ee es ee ee ee es ee ee ee ee ee 


READ THE CRITICAL VALUES FOR XG AND ZG FROM FILE DATA.DAT 


XGG(1)=0.0 

ZGG(1)=0.016358083 
DG te i=2197 

READ (21,*)XG,ZG 

XGG(I)=XG 

ZGG(I)=ZG 


DETERMINE [F] COEFFICIENTS 


D=( (M-YVDOT) * (IZZ-NRDOT) * (IXX-KPDOT) ) 
—( (M-YVDOT) * (IXZ-KRDOT) * (-IXZ-NPDOT) ) 
—( (M*XG-NVDOT) * (M*XG-YRDOT) * (IXX-KPDOT) ) 
+( (M*XG-NVDOT) * (IXZ-KRDOT) * (-M*ZG-YPDOT) ) 
+( (-M*ZG-KVDOT) * (M*XG-YRDOT) *(-IXZ-NPDOT) ) 
~ ((-M*ZG-KVDOT) * (IZZ-NRDOT) * (-M*ZG-YPDOT) ) 


A11=( (IZZ-NRDOT) * (IXX-KPDOT) ) - ( (IXZ-KRDOT) * (-IXZ-NPDOT) ) 
A12=( (-M*XG+YRDOT) * (IXX-KPDOT) ) + ( (IXZ-KRDOT) * (-M*ZG-YPDOT) ) 
A13=( (M*XG-YRDOT) *(-IXZ-NPDOT) ) - ( (IZZ-NRDOT) * (-M*ZG-YPDOT) ) 
A21=( (-M*XG+NVDOT) * (IXX-KPDOT) ) + ( (-M*ZG-KVDOT) * (-IXZ-NPDOT) ) 
A22=( (M-YVDOT) * (IXX-KPDOT) ) - ((-M*ZG-KVDOT) * (-M*ZG-YPDOT) ) 
A23=( (-M+YVDOT) * (-IXZ-NPDOT) ) + ( (M*XG-NVDOT) * (-M*ZG-YPDOT) ) 
A31=( (M*XG-NVDOT) * (IXZ-KRDOT) ) - ((-M*ZG-KVDOT) * (IZZ-NRDOT) ) 
A32=( (-M+YVDOT) * (IXZ-KRDOT) ) + ( (-M*ZG-KVDOT) * (M*XG-YRDOT) ) 
A33=( (M-YVDOT) * (IZZ-NRDOT) ) — ( (M*XG-NVDOT) * (M*XG-YRDOT) ) 


F(1,1)=(A11*YV*U+A12*NV*U+A13*KV*U) /D 
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12 
ia 


14 


18 


F(1,2)=(A11* (YR*U-M*U) +A12* (-M*XG*U+NR*U) +A13% (M*ZG*U+KR*U) ) /D 
F(1,3)=(A11*YP*U+A12*NP*U+A13*KP*U) /D 

F(1,4)=(A11* (W-B) +A12* (XG*W-XB*B) +A13* (-ZG*W+ZB*B) ) /D 
F(2,1)=(A21*YV*U+A22*NV*U+A23*KV*U) /D 

F(2,2)=(A21* (YR*U-M*U) +A22* (-M*XG*U+NR*U) +A23* (M*ZG*U+KR*U) ) /D 
F(2,3)=(A21*YP*U+A22*NP*U+A23*KP*U) /D 

F(2,4)=(A21* (W-B) +A22* (XG*W-XB*B) +A23* (-ZG*W+ZB*B) )/D 
F(3,1)=(A31*YV*U+A32*NV*U+A33*KV*U) /D 

F(3,2)=(A31* (YR*U-M*U) +A32* (-M*XG*U+NR*U) +A33* (M*ZG*U+KR*U) ) /D 
F(3,3)=(A31* YP*U+A32*NP*U+A33*KP*U) /D 

F(3,4)=(A31* (W-B) +A32* (XG*W-XB*B) +A33% (-ZG*W+ZB*B) )/D 
F(4,1)=0.0 

F(4,2)=0.0 

F(4,3)=1.0 

F(4,4)=0.0 


EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 


DO 11 K=1,4 
DO 12 J=1,4 
ASAVE(K, J) =F (K, J) 
CONTINUE 
CONTINUE 
CALL RG(4,4,F,WR,WI,1,YYY,IV1,FVi,IERR) 
WRITE(23,1007)WR(1) ,WR(2) , WR(3) , WR(4) 
CALL DSOMEGCIEV, WR, WI, OMEGA, CHECK) 
OMEGAO=OMEGA 
DO 5 J=1,4 
T(J,1)= YY¥ CO) TEV) 
T(J,2)=-YYY(J, IEV+1) 
CONTINUE 
IF (IEV.EQ.1.0) GO TO 13 
IF (IEV.EQ.2.0) GO TO 18 
IF (IEV.EQ.3.0) GO TO 14 
STOP 3004 
DO 6 J=1,4 
T(J,3)=YYY(J, 1) 
Day) =VYY (), 2) 
CONTINUE 
GO TO 17 
DO 19 J=1,4 
ees =Y Vy Cd, 1) 
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Q2 


al OQ 


©) 


8, 


13 


16 
at 


NO Ww 


T(J4) =a Gd 54) 
CONTINUE 
GO TO i7 
DO 16 J=1,4 
TJ, =r, 3) 
T(J ,4)=YYY (J, 4) 
CONTINUE 
CONTINUE 


NORMALIZATION OF THE CRITICAL EIGENVECTOR 
CALL NORMAL (T) 
INVERT TRANSFORMATION MATRIX 


DO 2 K=1,4 
DO 3 J=1,4 
TINV(K,J)=0.0 
TSAVE(K, J)=T(K, J) 
CONTINUE 
CONTINUE 
CALL DLUD(4,4, TSAVE,4,TLUD, IVLUD) 


DO 4 Jj=1,4 


TE SVE CI) EO. 0) Sree sco 
CONTINUE 
CALL DILU(4,4,TLUD, IVLUD,SVLUD) 
DO 8 K=1,4 

DO 9 J=1,4 

TINV(K, J)=TLUD (K, J) 

CONTINUE 

CONTINUE 


CHECK Inv(T) *A*T 

CALL MULT(TINV, ASAVE,T,A2) 
Pi=A2(3,3) 

P2=A2(4, 4) 

PEIGi=Pi 


PETG2=P2 
WRITE(25,1008)A2(3 51), A2(25 2) P12 
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©) 


DEFINITION OF Nij 


N1i1=TINV(1,1) 
N12=TINV(1,2) 
N13=TINV(1,3) 
N14=TINV (1,4) 
N21=TINV (2,1) 
N22-ciNy (22) 
N23=TINV (2,3) 
N24=TINV (2,4) 
Net-Tinves. 1) 
Ne2=TINV (3.2) 
N33=TINV (3,3) 
N34=TINV (3,4) 
N41=TINV(4,1) 
N42=TINV (4, 2) 
N43=TINV (4,3) 
N44=TINV (4,4) 


DEFINITION OF Mij 


M11=T(1,1) 
M12=T(1,2) 
M13=T (1,3) 
M14=T(1,4) 
M21=T(2,1) 
M22=T (2,2) 
M23=T (2,3) 
M24=T (2,4) 
M31=T (3,1) 
M32=T (3,2) 
M33=T (3,3) 
M34=1(3,4) 
M41=T(4,1) 
M42=T (4,2) 
M43=T (4,3) 
M44=T (4,4) 


DEFINITION OF Lij 
L1=YG* ((M31**2)+(M21**2) ) 


L2=IXY* (M31**2)+IYZ*M31*M21+YG*M21*M11 
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L3=IXY*M31*M21+IYZ* (M21**2) +M*YG*M11*M31+YG*W* (M41**2) 


& -YB*B* (M41**2) 


L4=YG* (2.0*M31*M32+2 .0*M21*M22) 
L5=2.0*IXY*M31*M32+1YZ* (M31*M22+M32*M21)+YG* (M21*M12+M22*M11) 
L6=IXY* (M31*M22+M32*M21)+2.0*1IYZ*M21#M22+M* YG* (M11*M32+M12*M31) 


& +2 .0* (YG*W-YB*B) *M41*M42 


& 


& 


& 


& 


& 


& 


& 


& 


& 
& 


L7=YG* ((M32**2)+(M22**2) ) 
L8=IXY* (M32**2)+IYZ*M32*M22+YG*M22*M12 
LO=IXY*M32*M22+1 YZ* (M22* *2) +M*YG*M12*M32+ (YG*W-YB*B) * (M42**2) 


L15=(A11/D)*L1+(A12/D) *L2-(A13/D)*L3 
L16=(A11/D)*L4+(A12/D) *L5-(A13/D) *L6 
L17=(A11/D) *L7+(A12/D) *L8-(A13/D) *L9 
L25=(A21/D)*L1+(A22/D) *L2-(A23/D)*L3 
L26=(A21/D) *L4+(A22/D) *L5-(A23/D) *L6 
L27=(A21/D)*L7+(A22/D) *L8-(A23/D) *L9 
L35=(A31/D) *L1+(A32/D) *L2-(A33/D) *L3 
L36=(A31/D) *L4+(A32/D) *L5- (A33/D) *L6 
L37=(A31/D)*L7+(A32/D) *L8- (A33/D) *L9 


C=CD/(6.0*GAMA) 

L1A=C* (EO* (M11**3)+3.O0*E1* (M11**2)*M21+3.0*E2*M11* (M21**2) 
+E3* (M21**3))+ (1.0/6.0) * (W-B) * (M41 **3) 

L2A=C* (E1* (M11**3) +3 .0*E2* (M11**2) *M21+3 .0*E3*M11* (M21**2) 
+E4* (M21**3) )+( (1.0/6.0) * (XG*W-XB*B) * (M41**3) 

L3A=(1.0/6.0)* (ZG*W-ZB*B) * (M41**3) 

L4A=C* (3. 0*EO* (M11**2) *M12+3 .0*E1* ((M11**2) #M22+2 .0*M11*M12*M21) 
+3 .0*E2* (M12* (M21**2) +2 .0*M21*M22*M11)+3.0*E3* (M21**2) *M22) 
+(1.0/6.0)* (W-B) *3.0* (M41**2) *M42 

L5A=C* (3.0*E1* (M11**2) *M12+3 .0*E2* ((M11**2) *M22+2 .0*M11*M12*M21) 
+3 .0*E3* (M12* (M21**2) +2 .0*M21*M22*M11)+3 .0*E4* (M21**2) *M22) 
+(1.0/6.0) * (XG*W-XB*B) *3.0* (M41**2) *M42 

L6A=(1.0/6.0)* (ZG*W-ZB*B) *3 .O* (M41**2) *M42 

L7A=C* (3.0*EO*#M11* (M12**2)+3 .0#E1* ((M12**2) *M214+2 .0*M11*M12*M22) 
+3 .0*E2* ((M22**2) *M11+2.0*M21*M22*M12)+3 .0*E3*M21* (M22**2) ) 
+(1.0/6.0)* (W-B) *3.0*M41* (M42**2) 

L8A=C* (3 .0*E1*M11* (M12**2)+3.0*E2* ((M12**2) *M21+2.0*M11*M12*M22) 
+3 .0*E3* ((M22**2) *M11+2 .0*M21*M22*M12)+3 .0*E4*M21* (M22**2) ) 
+(1.0/6.0) * (XG*W-XB*B) *3 .0*M41* (M42*+*«2) 

LOA=(1.0/6.0) * (ZG*W-ZB*B) *3 .0*M41* (M42**2) 

L10A=C* (EO* (M12**3) +3 .0*E1* (M11**2) *M21+3 .0*E2*M12* (M22**2) 
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& +E3* (M22**3) )+(1.0/6.0) * (W-B) * (M42**3) 

LiLA=C* (E1* (M12**3) +3 .O*E2* (M11**2) *M21+3 .0*E3*M12* (M22**2) 
& +E4* (M22**3) )+(1.0/6.0) * (XG*W-XB*B) * (M42**3) 
L12A=(1.0/6.0) * (ZG*W-ZB*B) * (M42**3) 


L11=(-A11/D) *L1A+(-A12/D) *L2A+(A13/D)*L3A 
L12=(-A11/D) *L4A+(-A12/D) *L5A+(A13/D) *L6A 
L13=(-A11/D) *L7A+(-A12/D) *L8A+(A13/D) *L9A 
Li4=(-A11/D) *L10A+(-A12/D) *L11A+(A13/D) *L12A 


L21=(-A21/D) *L1A+(-A22/D) *L2A+(A23/D) *L3A 
L22=(-A21/D) *L4A+(-A22/D) *L5A+(A23/D) *L6A 
L23=(-A21/D) *L7A+(-A22/D) *L8A+(A23/D) *L9A 
L24=(-A21/D) *L10A+(-A22/D) *L11A+(A23/D) *L12A 


L31=(-A31/D) *L1A+(-A32/D) *L2A+(A33/D) *L3A 
L32=(-A31/D) *L4A+(-A32/D) *L5A+(A33/D) *L6A 
L33=(-A31/D) *L7A+(-A32/D) *L8A+(A33/D) *L9OA 
L34=(-A31/D) *L10A+(-A32/D) *L11A+(A33/D) *L12A 


Rii=(N11*L11)+(N12*L21)+(N13*L31) 
Ri2=(N1i1*L12)+(N12*L22) +(N13*L32) 
R1i3=(N11*L13)+(N12*L23) +(N13*L33) 
R14=(N11*L14) +(N12*L24) +(N13*L34) 
R21=(N21*L11) + (N22*L21) +(N23*L31) 
R22=(N21*L12)+(N22*L22) + (N23*L32) 
R23=(N21*L13) + (N22*L23) +(N23*L33) 
R24=(N21*L14) + (N22*L24) +(N23*L34) 


P11=(N11*L15)+(N12*L25) +(N13*L35) 
P12=(N11*L16) +(N12*L26) +(N13*L36) 
P13=(N11*L17)+(N12*L27) +(N13*L37) 
P21=(N21*L15) +(N22*L25) +(N23*L35) 
P22=(N21*L16) + (N22*L26) +(N23*L36) 
P23=(N21*L17) +(N22*L27) +(N23*L37) 


EVALUATE DALPHA AND DOMEGA 
ZGR =ZGG(I) 
7Gle 766-1) 


ZG1 =ZGR 
XG1 =XGG(T) 
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~ CALL FMATRIX(ZG1,XG1,FF) 


CALL RG(4,4,FF,WR,WI,0,YYY,IV1,FV1,IERR) 
CALL DSTABL(DEOS,WR,WI, FREQ) 

ALPHR=DEOS 

OMEGR=FREQ 


ZGi =ZGL 
XG1 =XGG(I) 


CALL FMATRIX(ZG1,XG1,FF) 


CALL RG(4,4,FF,WR,WI,0,YYY,IV1,FV1,IERR) 
CALL DSTABL(DEOS, WR, WI, FREQ) 

ALPHL=DEOS 

OMEGL=FREQ 


DALPHA=(ALPHR-ALPHL) / (ZGR-ZGL) 
DOMEGA= (OMEGR-OMEGL) / (ZGR-ZGL) 


EVALUATION OF HOPF BIFURCATION COEFFICIENTS 


COEF1=(1.0/8.0)*(3.0*R11+R13+R22+3 .0*R24) 

COEF2=(1.0/8.0)*(3.0*R11+R23-R12-3.0*R14) 

WRITE (20,2001)XG,ZG,COEF1,DALPHA , OMEGAO, PEIG1 , PEIG2 
1 CONTINUE 


STOP 
2001 FORMAT (7E14.5) 
1007 FORMAT (4E14.5) 
1008 FORMAT (4E14.5) 
END 


oh 
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